Project Summary

This project aimed to look at the spatial variability of Symbiodinium clades C and D in the Kane’ohe Bay, O’ahu, Hawai’i population of Montipora capitata. We investigated the distribution of symbionts at scales ranging from location in the bay to location on an individual reef. We also looked at differences among reef types (fringing vs. patch), colony color morph (brown vs. orange) and depth. Heterogeneous mixtures of symbiont clades were considered in the analysis for spatial patterns. By investigating spatial variability of Symbiodinium, we furthered the understanding of stress-response potential in Kane’ohe Bay.

Library Packages

setwd("~/Symcap")
library(data.table)
library(devtools)
library(plyr)
library(reshape2)
library(popbio)
library(RgoogleMaps)
library(plotrix)
library(zoo)
library(rgdal)
library(car)
library(scales)

Import/Merge Field and qPCR Data

Coral_Data <- read.csv("Coral_Collection.csv")
Coral_Data$Depth..m. <- as.numeric(as.character(Coral_Data$Depth..m.))
source_url("https://raw.githubusercontent.com/jrcunning/steponeR/master/steponeR.R")
Mcap.plates <- list.files(path = "qPCR_data", pattern = "txt$", full.names = T)
Mcap <- steponeR(files = Mcap.plates, delim="\t",
                 target.ratios=c("C.D"), 
                 fluor.norm = list(C=2.26827, D=0), 
                 copy.number=list(C=33, D=3))
Mcap <- Mcap$result
Mcap <- Mcap[grep("+", Mcap$Sample.Name, fixed=T, invert = T), ]
Mcap <- Mcap[grep("NTC", Mcap$Sample.Name, fixed = T, invert = T), ]
Mcap <- Mcap[grep("PCT", Mcap$Sample.Name, fixed = T, invert = T), ]
colnames(Mcap)[which(colnames(Mcap)=="Sample.Name")] <- "Colony"
Mcap$fail <- ifelse(Mcap$C.reps < 2 & Mcap$D.reps < 2, TRUE, FALSE)
fails <- Mcap[Mcap$fail==TRUE, ]
Mcap <- Mcap[which(Mcap$fail==FALSE),]
Mcap$C.D[which(Mcap$C.reps<2)] <- -Inf
Mcap$C.D[which(Mcap$D.reps<2)] <- Inf
Mcap <- Mcap[with(Mcap, order(Colony)), ]
Mcap$propC <- Mcap$C.D / (Mcap$C.D+1)
Mcap$propD <- 1-Mcap$propC
Mcap$propD[which(Mcap$C.D==-Inf)] <-1
Mcap$propC[which(Mcap$C.D==-Inf)] <-0
Mcap$propD[which(Mcap$C.D==Inf)] <-0
Mcap$propC[which(Mcap$C.D==Inf)] <-1
Mcap$Dom <- ifelse(Mcap$propC>Mcap$propD, "C", "D")
Symcap<-merge(Coral_Data, Mcap, by="Colony", all=T)
Symcap <- Symcap[with(Symcap, order(Colony)), ]
Symcap$Mix <- factor(ifelse(Symcap$propC>Symcap$propD, ifelse(Symcap$propD!=0, "CD", "C"), ifelse(Symcap$propD>Symcap$propC, ifelse(Symcap$propC!=0, "DC", "D"), NA)), levels = c("C", "CD", "DC", "D"))
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")

Adjust Depth by Mean Sea Level

To account for the influence of tides, depth was adjusted according to the difference in sea level from the mean sea level on each collection day at 6-minute intervals. Mean sea level was obtained from NOAA tide tables for Moku o Lo’e.

JuneTide=read.csv("Station_1612480_tide_ht_20160601-20160630.csv")
JulyTide=read.csv("Station_1612480_tide_ht_20160701-20160731.csv")
Tide<-rbind(JuneTide, JulyTide)
Tide$Time <- as.POSIXct(Tide$TimeUTC, format="%Y-%m-%d %H:%M:%S", tz="UTC")
attributes(Tide$Time)$tzone <- "Pacific/Honolulu"
Symcap$Time2 <- as.POSIXct(paste(as.character(Symcap$Date), as.character(Symcap$Time)), format="%m/%d/%y %H:%M", tz="Pacific/Honolulu")
Symcap$Time=Symcap$Time2
Round6 <- function (times)  {
  x <- as.POSIXlt(times)
  mins <- x$min
  mult <- mins %/% 6
  remain <- mins %% 6
  if(remain > 3L) {
    mult <- mult + 1
  } else {
    x$min <- 6 * mult
  }
  x <- as.POSIXct(x)
  x
}
Symcap$Time.r <- Round6(Symcap$Time)
Tide$Time.r <- Tide$Time
merged<-merge(Symcap, Tide, by="Time.r", all.x=T)
merged$newDepth <- merged$Depth..m.- merged$TideHT

Dominant Symbiont by Reef Area

Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Dom, Symcap$Reef.Area)
chisq.test(results)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 136.26, df = 1, p-value < 2.2e-16
prop.table(results, margin = 2)
##    
##         Slope       Top
##   C 0.7767857 0.3294574
##   D 0.2232143 0.6705426
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Symbmiont Proportion")
legend("topright", legend=c("C", "D"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)

Dominant Symbiont by Color Morph

results=table(Symcap$Dom, Symcap$Color.Morph)
chisq.test(results)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 164.96, df = 1, p-value < 2.2e-16
prop.table(results, margin = 2)
##    
##         Brown    Orange
##   C 0.8896321 0.4103194
##   D 0.1103679 0.5896806
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Color Morph", ylab = "Symbiont Proportion")
legend("topright", legend=c("C", "D"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)

Symbiont Community Composition by Color Morph

results=table(Symcap$Mix, Symcap$Color.Morph)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 167.44, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##     
##            Brown      Orange
##   C  0.762541806 0.361179361
##   CD 0.127090301 0.049140049
##   DC 0.107023411 0.570024570
##   D  0.003344482 0.019656020
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Color Morph", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)

Symbiont Community Composition by Reef Area

Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Mix, Symcap$Reef.Area)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 138.97, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##     
##            Slope         Top
##   C  0.678571429 0.275193798
##   CD 0.098214286 0.054263566
##   DC 0.214285714 0.651162791
##   D  0.008928571 0.019379845
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Reef Area", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)

Symbmiont Community Composition by Dominant Symbiont

results=table(Symcap$Mix, Symcap$Dom)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 707, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##     
##               C          D
##   C  0.86635945 0.00000000
##   CD 0.13364055 0.00000000
##   DC 0.00000000 0.96703297
##   D  0.00000000 0.03296703
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray 10", "gray 85", "gray 40", "gray100"), xlab = "Dominant Symbiont", ylab = "Symbiont Mixture Proportion")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray85", "gray40", "gray100"), inset = c(-.2, 0), xpd = NA)

Proportion of D When Present in Mixture

propD <- merged$propD[which(merged$propD > 0 & merged$propD < 1)]
hist(propD, xlab = "Proportion of Clade D", ylab = "Number of Samples", main = "", col = "gray75")

Color Morph by Reef Area

Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Color.Morph, Symcap$Reef.Area)
chisq.test(results)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 81.109, df = 1, p-value < 2.2e-16
prop.table(results, margin = 2)
##         
##              Slope       Top
##   Brown  0.5523385 0.2015504
##   Orange 0.4476615 0.7984496
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Color Morph Proportion")
legend("topright", legend=c("Brown", "Orange"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)

Color Morph by Reef ID

results=table(Symcap$Color.Morph, Symcap$Reef.ID)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 57.007, df = 25, p-value = 0.0002667
prop.table(results, margin = 2)
##         
##                 10        13        14        18        19        20
##   Brown  0.4000000 0.4000000 0.5333333 0.2333333 0.3666667 0.6285714
##   Orange 0.6000000 0.6000000 0.4666667 0.7666667 0.6333333 0.3714286
##         
##                 21        25        26        30        37        38
##   Brown  0.4333333 0.6000000 0.4242424 0.1333333 0.8000000 0.2000000
##   Orange 0.5666667 0.4000000 0.5757576 0.8666667 0.2000000 0.8000000
##         
##                 42        46         5      Deep     F1-46     F2-25
##   Brown  0.5142857 0.5000000 0.3666667 0.5172414 0.4000000 0.6000000
##   Orange 0.4857143 0.5000000 0.6333333 0.4827586 0.6000000 0.4000000
##         
##              F3-18     F4-34     F5-34 F6-Lilipuna     F7-18     F8-10
##   Brown  0.2500000 0.4400000 0.2500000   0.5000000 0.2000000 0.4000000
##   Orange 0.7500000 0.5600000 0.7500000   0.5000000 0.8000000 0.6000000
##         
##               F9-5      HIMB
##   Brown  0.4000000 0.6000000
##   Orange 0.6000000 0.4000000

Dominant Symbiont by Depth

merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
Dom1 <- subset(merged, !is.na(newDepth) & !is.na(Dominant))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
logi.hist.plot(Dom1$newDepth, Dom1$Dominant, boxp = FALSE, type = "hist", col="gray", xlabel = "Depth (m)", ylabel = "", ylabel2 = "")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "C                             D", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of clade C Symbiont", line = 3, cex = 1)

merged$DepthInt <- cut(merged$Depth..m., breaks = 0:13)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.2, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Dominant~Depth..m., family = "binomial", data = merged)
fitted <- predict(results, newdata = list(Depth..m.=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Dominant Symbiont Proportion")

Bars indicate relative frequency of clade C vs. D dominance at 1m depth intervals when pooling all samples collected.

Symbiont Community Composition by Depth

merged$Mixture <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 1, 0)
merged$Mixture2 <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 0, 1)
results=glm(Mixture~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Mixture
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                       589     817.24              
## newDepth  1   75.664       588     741.57 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
results=table(merged$Mixture2, merged$DepthInt)
results
##    
##     (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
##   0   172    43    30    34    18    11     3     2     5      2       1
##   1    66    41    60    99    46    20    10    17    12      6       6
##    
##     (11,12] (12,13]
##   0       1       0
##   1       0       1
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Mix", "Non-Mix"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.23, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Mixture~Depth..m., family = "binomial", data = merged)
fitted <- predict(results, newdata = list(Depth..m.=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Mixture Proportion")

Color Morph by Depth

merged$Color <- ifelse(merged$Color.Morph=="Orange", 1, 0)
results=glm(Color~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
Color <- subset(merged, !is.na(newDepth) & !is.na(Color))
logi.hist.plot(independ = Color$newDepth, depend = Color$Color, type = "hist", boxp = FALSE, ylabel = "", col="gray", ylabel2 = "", xlabel = "Depth (m)")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "Brown                       Orange", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)

merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(-.22, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~Depth..m., family = "binomial", data = merged)
fitted <- predict(results, newdata = list(Depth..m.=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Color Morph Proportion")

Bars indicate relative frequency of Brown vs. Orange color morph dominance at 1m depth intervals when pooling all samples collected.

Dominant Symbiont by Color Morph and Reef Area

merged$Reef.Area <- ifelse(merged$Reef.Area!="Top", yes = "Slope", no = "Top")
table(merged$Dom, merged$Color.Morph, merged$Reef.Area)
## , ,  = Slope
## 
##    
##     Brown Orange
##   C   226    122
##   D    21     79
## 
## , ,  = Top
## 
##    
##     Brown Orange
##   C    40     45
##   D    12    161
model=aov(Dominant~Color.Morph*Reef.Area, data = merged)
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: Dominant
##                        Sum Sq  Df F value    Pr(>F)    
## Color.Morph            21.329   1 134.208 < 2.2e-16 ***
## Reef.Area              14.489   1  91.168 < 2.2e-16 ***
## Color.Morph:Reef.Area   1.780   1  11.201 0.0008612 ***
## Residuals             111.565 702                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dominant Symbiont by Color Morph and Depth

Because and interactive effect of reef area and color morph was observed and slope is indicative of a depth gradient, the interaction between depth and color morph was tested.

model1=lm(Dominant~Color.Morph*newDepth, data = merged)
Anova(model1, type = 2)
## Anova Table (Type II tests)
## 
## Response: Dominant
##                      Sum Sq  Df F value    Pr(>F)    
## Color.Morph          24.863   1 154.391 < 2.2e-16 ***
## newDepth             12.236   1  75.981 < 2.2e-16 ***
## Color.Morph:newDepth  2.592   1  16.097 6.798e-05 ***
## Residuals            94.369 586                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dominant Symbiont by Depth and Reef Type

model3=aov(Dominant2~Depth..m.*Reef.Type, data = merged)
Anova(model3, type = 2)
## Anova Table (Type II tests)
## 
## Response: Dominant2
##                      Sum Sq  Df  F value    Pr(>F)    
## Depth..m.            24.214   1 120.3238 < 2.2e-16 ***
## Reef.Type             0.052   1   0.2569  0.612431    
## Depth..m.:Reef.Type   1.727   1   8.5815  0.003506 ** 
## Residuals           141.272 702                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df <- subset(merged, Reef.Type=="Patch")
results=glm(Dominant2~Depth..m., family = "binomial", data = df)
newdata <- list(Depth..m.=seq(0,12,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,12,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
results=glm(Dominant2~Depth..m., family = "binomial", data = df)
newdata <- list(Depth..m.=seq(0,12,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,12,0.01), col="brown1", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Clade C Symbiont", line = 3, cex = 1)
legend("topright", legend=c("Patch", "Fringe"), fill=c("dodgerblue3", "brown1"), inset = c(-.22, 0), xpd = NA)

Depth at which Orange switches from D to C Dominance

While brown was always C-dominated, orange was observed to switch from D to C dominance. The depth threshold at which this switch occurs was analyzed here.

threshdepth <- function(color) {
  df <- subset(merged, Color.Morph==color)
  plot(df$Dominant2~df$newDepth, xlab="Depth (m)", ylab = "Probability of Clade C Symbiont",
       main=color)
  abline(h = 0.5, lty=2)
  results=glm(Dominant2~newDepth, family = "binomial", data = df)
  pval <- data.frame(coef(summary(results)))$`Pr...z..`[2]
  mtext(side=3, text=pval)
  newdata <- list(newDepth=seq(0,12,0.01))
  fitted <- predict(results, newdata = newdata, type = "response")
  lines(fitted ~ seq(0,12,0.01))
  thresh <- ifelse(pval < 0.05,
                   newdata$newDepth[which(diff(sign(fitted - 0.5))!=0)], NA)
  return(thresh)
}
sapply(levels(merged$Color.Morph), FUN=threshdepth)

##  Brown Orange 
##     NA    2.7
df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,12,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,12,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,12,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,12,0.01), col="sienna", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Clade C Symbiont", line = 3, cex = 1)
legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(-.22, 0), xpd = NA)

Color Morph by Reef Type and Depth

model2=aov(Color~newDepth*Reef.Type, data = merged)
Anova(model2, type = 2)
## Anova Table (Type II tests)
## 
## Response: Color
##                     Sum Sq  Df F value    Pr(>F)    
## newDepth             6.180   1 27.6214 2.069e-07 ***
## Reef.Type            0.132   1  0.5908  0.442436    
## newDepth:Reef.Type   1.651   1  7.3813  0.006785 ** 
## Residuals          131.109 586                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df <- subset(merged, Reef.Type=="Patch")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,12,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,12,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,12,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,12,0.01), col="brown1", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)
legend("topright", legend=c("Patch", "Fringe"), fill=c("dodgerblue3", "brown1"), inset = c(-.22, 0), xpd = NA)

Figures

Coral Collection Map

KB <- c(21.46087401, -157.809907) 
KBMap <- GetMap(center = KB, zoom = 13, maptype = "satellite", SCALE = 2, GRAYSCALE = FALSE)
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
XY <- subset(XY, Reef.ID!="37")
par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude, col=153, pch=21, bg="#FF7F50", lwd=2)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)
par(new=T, mar=c(11.8,19,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83") 
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83")) 
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()

Depth Influence on Dominant Symbiont and Color Morph

par(mfrow=c(3,1))

merged$DepthInt <- cut(merged$Depth..m., breaks = 0:13)
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(2, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(2.1, 4, 2, 6))
results=glm(Dominant~Depth..m., family = "binomial", data = merged)
fitted <- predict(results, newdata = list(Depth..m.=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="Probability of D-Dominance")

merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(3, 4, 1, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)), 
        xlab = "", ylab = "Probability of Orange-Dominance",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(3.1, 4, 1, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~Depth..m., family = "binomial", data = merged)
fitted <- predict(results, newdata = list(Depth..m.=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="")

df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 0, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="Depth (m)", ylab="Probabilty of D-Dominance", axisnames=FALSE, xaxs = "i", yaxs = "i")
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0, 11, 0.01), col="sienna", lwd=3)
#legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(-.13, 0), xpd = NA)

results=table(merged$Dominant2, merged$Color.Morph)
chisq.test(results)
prop.table(results, margin = 2)
par(new=T, mar=c(10, 10, .5, 6.3))
barplot(prop.table(results, margin = 2), col = c(alpha("red", 0.35), alpha("blue", 0.35)), xlab = "", ylab = "", yaxt = 'n')

#legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.35), alpha("red", 0.35)), inset = c(0, 0), xpd = NA)